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The solutions of the one-dimensional homogeneous nonlinear Boltzmann equation are studied in 
the QE-limit (Quasi-Elastic; infinitesimal dissipation) by a combination of analytical and numerical 
techniques. Their behavior at large velocities differs qualitatively from that for higher dimensional 
systems. In our generic model, a dissipative fluid is maintained in a non-equilibrium steady state 
by a stochastic or deterministic driving force. The velocity distribution for stochastic driving is 
regular and for infinitesimal dissipation, has a stretched exponential tail, with an unusual stretching 
exponent bqE = 2b, twice as large as the standard one for the corresponding d-dimensional system 
at finite dissipation. For deterministic driving the behavior is more subtle and displays singular- 
ities, such as multi-peaked velocity distribution functions. We classify the corresponding velocity 
distributions according to the nature and scaling behavior of such singularities. 
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I. INTRODUCTION 



A. Background and outline 



The model of inelastic hard spheres is one of the most simple frameworks to describe granular gases (see e.g. 
0, [H, [H for reviews and further references). The contraction of phase space due to dissipative collisions leads to a non- 
equilibrium behavior that is markedly different from that of equilibrium systems (non-Gaussian velocity distributions, 
counter- intuitive hydrodynamics, breakdown of kinetic energy equipartition etc [H, 0, HI ) • In this paper, we study in 



detail the limit of quasi-elasticity 
been the subject of some interest 



with particular emphasis on one-dimensional systems, that have already 

The kinetic description is provided by the nonlinear Boltzmann equation. As we are interested in the velocity 
statistics of dissipative gases, we will restrict ourselves to homogeneous and isotropic solutions. Any spatial dependence 
will therefore be discarded. The time evolution of the velocity distribution function F(v,t) is then governed by the 
following eauation|13. lD|, 



d t F(v)+TF = I{v\F) = 



d-wg" 



-^ I F(v**)F(w**) - F(v)F(w) 



(1) 



Here TF represents the action of a driving mechanism, that injects energy into the system, and counterbalances the 
energy dissipated by inelastic collisions. Consequently, the system is expected to reach a non-equilibrium steady- 
state. In the equation above g = v — w denotes the relative velocity of colliding particles with g = |g|, J (• ■•) = 
(l/fid) J dn(- ■ ■ ) is an angular average over the surface area = 2Tr d ' 2 /T(^d) of a d-dimensional unit sphere, and 
g v models the collision frequency. Note that in one dimension the integral f n is absent. We have absorbed constant 
factors in the time scale. Here (v**, w**) denote the restituting velocities that yield (v, w) as post-collisional velocities, 
i.e. 



= v-i(l 



"(/Kg- 11 ) 11 



= w + ±(l + a 1 )(g-n)n 



(2) 



where (unit) vector n is parallel to the line of centers of the colliding particles. Note that nn is replaced by 1 in one 
dimension. The direct collision law is obtained from (|2|) by interchanging pre- and post-collision velocities and by 
replacing ao - * l/ a o where ao < 1 is the restitution coefficient. Each collision leads to an energy loss proportional to 
(1 — <Xq). Elastic collisions therefore correspond to ao = 1. 
In this article, we will consider the source term in |T]) to be of the form 



TF = d ■ (aF) - Dd 2 F = jd ■ (vv e F) ~ Dd 2 F, 



(3) 



where a = vv e is a negative friction force, d = d/dv is the gradient in v-space, and 7 and D are positive constants. 
Two situations will be addressed : (7 = 0, D > 0) or (7 > 0, D — 0). They correspond respectively to stochastic 



2 



White Noise (WN), or to deterministic nonlinear Negative Friction (NF). While the WN driving mechanism has been 
extensively studied @, [l5|, [l6| , the Negative Friction has been introduced more recently [l3| . The continuous exponent 
9 > selectively controls the energy injection mechanism. Schematically, increasing the value of 9 corresponds to 
injecting more ener gy i n the large velocity tail of the distribution. However two special values have been studied in 
the past [l2l Il3l. Hi] lid. fl7l fl^ . i.e. (i) the Gaussian thermostat (9 = 1), which is equivalent to the homogeneous free 
cooling state, where the system is unforced and the possibility of spatial heterogeneities discarded (see e.g. [Ill), anc ^ 
(ii) the case 9 = 0, referred to as 'gravity thermostat' [16], or as 'negative solid friction' [H, [3| . 

We emphasize that 7 and D are irrelevant constants that can be eliminated (see below), whereas the exponents 9 
and v are fundamental quantities for our purposes. As it appears in Eq. dJ), v governs the collision frequency of the 
system: v = 1 corresponds to hard-sphere like dynamics and v = to the so-called Maxwell model (I2I [13. 120LT2TL l22| . 
Our collision kernel generalizes these two cases to a general class of repulsive power law potentials, where v is related 
to the power law exponent and the dimensionality (see [Hj]). 

Under the action of the driving term T , the solution of (JTJ evolves towards a non-equilibrium steady state. We 
will be interested in the properties of the corresponding velocity distribution F(v), in the limit where ao — ► 1 — . 
Since the limit ao — ► 1 turns out to be singular in one dimension, attention must be paid to the fact that the value 
ao = 1 (elastic interactions) has to be analyzed separately. Indeed, when ao = 1 in one dimension, the collision law 
(|2|) simply corresponds to an exchange of particle labels. So the initial velocity distribution does not evolve in time. 
On the other hand, the steady state velocity distribution at any ao < 1 is independent of the initial condition. In 
dimensions higher than 1, this property holds for all values of ao < 1. In other words, the one dimensional situation 
with ao < 1 exhibits universal features, unlike its elastic counterpart. The quasi-elastic limit ao — > 1~ is therefore 
peculiar since a point with no universal properties (ao = 1) is approached via a 'universal' route (ao < 1). The 
resulting behavior of F shows some surprising features, that may be considered as mathematical curiosities, but are 
analytically challenging. It also turns out that they are numerically difficult to study. The numerical study relies on 
the DSMC (Direct Simulation by Monte Carlo) algorithm [23| , which allows us to obtain an exact numerical solution 
of the Boltzmann equation. As ao approaches 1, the memory of the initial conditions lasts for longer and longer times 
so that the computer time needed to reach the non-equilibrium steady state increases and simulations become more 
and more time-consuming. 

The behavior of our system is somewhat simpler in the case where energy is injected by a stochastic force (WN), 
and we start by analyzing this driving mechanism in section [Til It will be shown that the regular high energy tail of 
F(y) ~ exp[— v b ], which holds in any space dimension [12. Il3l Il4j. is preempted by a quasi-elastic tail characterized 
by a different stretching exponent 6qb = 2b. This is a signature of the non-commutativity of the limits: v — > 00 and 
ao — > l - . The case of driving through negative friction will be addressed in detail in section HTT1 As already observed 
for the homogeneous cooling state of inelastic hard rods [3, 0, [13, [H| , the velocity distribution becomes singular in 
the QE- limit, where it may approach a multi-peaked solution, and not a Gaussian. Starting from a small-inelasticity 
expansion for the collision operator I(v\F) in (TTJ, we characterize the scaling behavior. By a combination of analytical 
work and numerical evidence, we propose a classification of the different types of limiting velocity distributions, several 
of which correspond to new types of solutions of the nonlinear Boltzmann equation. 

B. Preliminary remarks 

We start by introducing some notations and summarizing a few results fl3l j that are relevant for our study. In the 
subsequent analysis, it is convenient to introduce the variables p = (1 + ao)/2, q = (1 — ao)/2 so that p + q = 1, 
and to measure the velocities in units of the r.m.s. velocity. We study steady states and introduce a rescaled velocity 
distribution /(c) such that F(v) = Vq f(v/vo) where v$ is the r.m.s velocity, c = v/vo, and d is the number of 
spatial dimensions. By definition, J dcf = 1 and the normalization chosen reads J dec 2 f = d/2. After inserting 
the scaling ansatz in Eq.([T]), we have shown in [l3| that a stable steady state for WN driving is reached provided 
bwN = 1 + v/2 > 0, and for NF driving provided &jv_f = v + 1 — 9 > 0. When b < 0, the non-equilibrium steady state 
is unstable, i.e. it is a repelling fixed point of the dynamics. Our analysis should therefore be limited to the cases 
v > -2 for WN and to v > 9 - 1 for NF. 

We have shown in [13j ] that the quantity b introduced above not only separates stable from unstable situations, but 
also governs the high energy tail of f{c). In marginal cases where b vanishes, / has a power-law tail. The freely cooling 
(9 = 1) Maxwell model [y = 0) provides an illustration that has been discussed in (20I. liHl. I22]] . On the other hand, 
when b > 0, / has a stretched exponential tail so that In /(c) oc — c b at large c. This result holds in any dimension. In 
d = 1, the correspon ding tail may however be 'masked' when ao is close to unity, a phenomenon already observed for 
hard rods (y = 1) in [17]]: the behavior In /(c) oc — c b holds for c > c*(ao), where c*(ao) is an ao-dependent threshold. 
In the QE-limit, the threshold value c*(ao) ^00. As a consequence, considering the limit of large c at any finite ao, 
the standard behavior with exponent b is observed. Alternatively, taking first the limit ao — > 1~ at fixed c , new tails 
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may appear. It is the purpose of the present paper to study their properties. In addition, in the case of NF driving, 
/(c) becomes singular in the quasi-elastic limit. Our goal will then be to understand the underlying scaling behavior 
and to propose a classification of the various limiting shapes for the velocity distribution. 

II. WHITE NOISE (WN) DRIVING 

Unless explicitly stated, we limit ourselves to d = 1. In the case of WN driving the integral equation |T]) for the 
scaling form, /(c) = vqF(cvq), becomes with the help of the relation I(v\F) — I(c\f): 

I(c\f) = -Dv-'^fic) = -± pqKv f»(c). (4) 

The second equality has been obtained by applying J dcc 2 (- ■ • ) to the first equality, see (I3|), yielding 

Dv^- 1 = -I« c 2 /(c|/))) = I W «|c- Ci r 2 )) = \ PQ k v , (5) 

where the double brackets denote an average with weight /(c)/(ci). Similarly, simple brackets denote an average with 
weight /(c). 

It is next convenient to treat the 1-D collision term, 

J(c|/)= [d Cl \c- Cl \» [a^- 7(c")/(cD-/(c)/(ci)] , (6) 



by using the inverse transformation, c = qu + pc\* and c\ = pu + 9c**, where u = c** and \c — c\\ = ao\c — u\/p. We 
obtain 



J(c|/)= / du\c-u\ v 



p-»- l f{u)f(^—^-yf(u)f(c)y (?) 

In the quasi-elastic limit q — > + , p — > 1~, we perform the small-q expansion, following Refs.[3, H, @], 

-^r/ ( £ ^ ?H ) = (1 + + + <?( c - «)/'( c ) + <^ 2 )- (8) 

Then we find to 0{q) included, 

J( C |/) = q /" d«/(«)|c-t»|"[(c-ti)/'(c) + (l/ + l)/(c)] 

- qfj(c) J du\c~u\»{c-u)f{u). (9) 

Note that these results hold irrespective of the driving mechanism. The second equality can be verified by evaluating 
the derivative. Inserting of ^ into @ allows us to integrate ((H) once, yielding: 



/'(c) + /(c)(2//c) J du\c - u\ v {c - u)f(u) = 0. (10) 
This equation can be integrated once more to obtain the implicit equation, 

/(c) = Cexpf- 2 / du\c-u\^f(u) 

-Cexpj^^-J (coo), (11) 

where C is an integration constant. The large velocity tail immediately follows. We recover a result already obtained 
in 0, QjJ for hard rods [y = 1), as well as one for Maxwell models {y = 0) in [l8j], and we note that in the QE- limit 
the exponent 6qs = 2bwN = v + 2. On the other hand, for any finite e = 1 — ao, the large c-behavior is given by 
In /(c) cx — c b with b = 1 + v/2 [T3|. These two facts show that the limits c — > 00 and e — ► do not commute, 

/(e, c) ~ exp(— A c 1+!y / 2 ) (standard tail: large c, fixed e 7^ 0) 

/(£,c)~ exp^Ac^ 2 ), (QE-tail: small e, fixed c). (12) 
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or, in a more rigorous formulation 



In f(e,c) 

lim, lira — — — 'j— = C (standard tail) 

lim !™ + = C ' (QE-tail). (13) 

c— too e— >Cr c ^ 



We have successfully tested these predictions against Direct Simulation Monte Carlo data [23|. Figure Q] shows that 
the standard tail with exponent &qb = 1 + v/2 is clearly observed for ao = (see inset), while the QE-tail with 
exponent v + 2 applies for ao = 0.995 (see main frame). Figure [2] conveys a similar message, and shows further that 
in two dimensions, Gaussian behavior is recovered, as expected, when ao — > 1 (see the inset of the right hand side 
figure). This illustrates the qualitatively different nature of the quasi-elastic limit in d = 1, as opposed to higher 
dimensions. 




c 



FIG. 1: Rescaled velocity distribution /(c) obtained by solving numerically the 1-D Boltzmann equation |T]) with WN driving 
by the DSMC technique [231 ] for the strongly interacting so-called 'very hard particles' model (v — 2) 19j). Two extreme cases 
have been simulated, one close to perfect elasticity (ao = 0.995) where the QE-tail In / vs. c" +2 (main frame) is visible, and 
the other at complete inelasticity (qq = 0), where the standard tail In / vs. c 1+v ^ 2 (inset) is visible. 




5 



III. NEGATIVE FRICTION (NF) DRIVING 

The scaling equation for the analog of ([4| with NF driving becomes, 

™^sW w 'l^i (+r,,) ' (14) 

Here I(c\f) takes the form ©, and the analog of ([5]) has been used to eliminate 7. Moreover, fig+i = (|c| e+1 ) and 
k v = ((|c— ci\ v+2 )). We therefore have to 0(q) included, 

f(c) [ dc x \c-c x \»{c-c l )f{c l ) = -^c\c\ e - l f{c). (15) 
J 2^0+1 

Previous studies of the free cooling regime of hard rods (9 = l,z/ = 1) have shown that /(c) becomes singular 
when c<o — * 1 ~, and evolves into two symmetric Dirac peaks [3, [^, [l?]]- It is indeed easy to check that /(c) = 
[S(c + a) + 8(c — a)]/2 is a solution of (|15[) h with a = 1/V2 as required by our choice of normalization, (c 2 ) = 1/2, 
and K„/(2/z e+ i) = 2"a l/+1 ~ e . Figure [3] shows that the approach to such a solution may be observed in the numerical 
simulations for (0, v) 7^ (1, 1). In addition, one can check that /(c) = A[S(c + a) + 5(c — a)] + B6(c) is also a solution 
of Eq. (fT5|) . provided that 2A + B = 1 and AAa 2 = 1 to enforce normalization. The DSMC results may indeed display 
such a three-peak structure (see Figure|3]). Note that the parameters corresponding to Figures [3] and 0] are quite close: 
= (1.0,1.2) and (1.1, 1.3) respectively. 
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FIG. 3: Velocity distribution obtained from DSMC simulations with NF driving for (6», v) = (1.0, 1.2) in 1-D. As e = 1-ao -> + , 
/(c) becomes increasingly peaked around c = ±c* = ±l/v2, and ultimately evolves into two symmetric Dirac distributions. 
All distributions displayed here and in other figures are such that J dcf — 1, f dccf = and J dec 2 f = 1/2. 

However, upon changing the parameters 9 and v, it appears that more complex shapes can be observed: /(c) may 
evolve towards a 4-, 5-, 6-peaks form, or other structures such as displayed in Figure [5] where /(c) seems to diverge 
at some points when «o — ► l - , with nevertheless a finite support. The diversity of the various velocity distributions 
obtained numerically calls for a rationalizing study. By a combination of analytical work and numerical evidence, we 
will propose below a classification of the different possible limiting velocity distributions. In addition, in the double 
peak case, two natural questions will be addressed: Are the peaks exemplified in Figure [3] self-similar ? If so, what is 
their shape ? 




A. The double peak scenario : structure and scaling 



We start by looking for scaling solutions to Eq. |T]), and restrict our analysis to the limiting form with the symmetric 
double peak. As in [TtJ we take /(c) of the doubly peaked form, 



^ = 4£ 



l + bc 
2E 



l-6c 
2E 



(16) 



where the width E = is expected to vanish when q — > + and 6, u, together with the function ip are unknowns. 
We impose (1)^ = ftp = 1, and we choose (x)^ = 0, which together with the condition (c 2 ) = 1/2, implies 



6 



10 



jj — e = 10~ 2 




: 


ji -- 1- to ; 






|; ■■■ t- to ; 








I 
ii 

H 

fj 

r.\ 


Jr 




-^Pf- — 





-1 -0.75 -0.5 -0.25 0.25 0.5 0.75 1 

c 

FIG. 4: Same as Figure[3]for slightly different parameters (9, v) = (1.1, 1.3). The velocity distribution now reaches a three-peak 
structure when e = 1 — ao — > + . 
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FIG. 5: Velocity distribution obtained from DSMC for (8, v) — (1.0,0.5). These parameters correspond to the border of the 
"Zoo" region shown in Figure [6] 



6 2 = 2(1 + 4E 2 (x 2 }^,), where normalization requires that b — > (V%)~ as q — > + . We note that asymmetric forms with 
ip(x) ^ ip{—x) may be realized (see [I?} and later sections). The ansatz (|T6)) allows us to resolve the structure of the 
Dirac peaks, shown in Figure [31 and to identify the type of self-similar behavior involved. However to this end, we 
need to expand the collision operator I(c\f) to second order in the inelasticity e — 1 — ao = 2g. Restriction to first 
order, as done in enables us to show that the double Dirac form is a solution of the Boltzmann equation, but 
does not allow us to impose constraints on the shape of ip and on the scaling exponent lu. Technical details can be 
found in the appendix, where it is shown that the equation fulfilled by ij}(x) reads, 

+ 2E 2 {v + 1 - 29)xip(x) + 2E v+2 tj){x){\x - y\ u {x - y))^ 
-E z v(y + l)^(x)(x 2 + (y 2 )^) + 2E 3 [{y + l)(y + 2) - 29(6 + 1)] (y 2 )^ ^(x) 

+4E 3 6(9 - l)x 2 ip{x) + 2E» +3 i;{x)((\y - z\ u+2 )) + 0{EqiP) = 0, (17) 

where y, z are dummy variables, and terms of O(Eqip) have been neglected. This relation involves terms of various 
orders in inelasticity. Given that E — q u , these terms are 0(q a ) with a = 1, 2u>, (2 + v)uu, 3a;, (3+v)w,<j + 1. Depending 
on the values of the parameters 9 and v one has to distinguish various possibilities to characterize the phase diagram, 
i.e. the physically allowed region of the (9, ^-parameter space. The stability criteria for the steady state (see [l3j]) 
constrain the phase diagram in Figure [S] to be inside the region [0 > 0, v > 9 — 1]. 

1. Case v > (a- and /3-scalings). 

Type a-scaling: The terms of order q and q 2u) are the dominant ones in Eq. (jTTJ) . So, u> = 1/2 and 

ij>'(x) = -(v+l-26)xip(x). (18) 
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FIG. 6: Phase diagram in the quasi-elastic limit. Each dot corresponds to a set of numerical simulations at smaller and 
smaller values of e = 1 — ao, performed to check the validity of the scaling behaviors summarized in Table [I] The stability 
of the non-equilibrium steady state requires that v > 9 — 1, while on the diagonal v — 9 — 1 the velocity distribution has a 
power-law tail fl3fl . Here a-, a'- and /3-scalings are associated with a distribution with two Dirac peaks. For 9 > 1 we observed 
numerically a solution with three Dirac peaks, while there does not seem to be a simple common feature for the distributions 
in the triangular "Zoo" region. An example of type a-scaling is given in Figure [7] (See also Figure [8] for type f3, and Figure [9] 
for type a'). Figures [12] and [13] together with Figure [5] above give an overview of several scaling shapes encountered in the Zoo 
region. The cross at (1,1) corresponds to the homogeneous cooling state of inelastic hard rods, explored in [l7j . 




x 

FIG. 7: Scaling behavior of type a. Plots of e w /(c) vs. x = (jcj — c*)/e^ at 9 = 1 for various inelasticities e — 1 — ao with 
io = 1/2 as predicted. Here c* = l/v2 corresponds to the peak of the distribution. Open symbols correspond to v — 1.5 and 
filled symbols are for v = 1.2 (same as in Figure EJ). The inset shows the same results on a linear-log scale. In each case, the 
prediction of Eq. (|19[1 f° r the scaling function is shown by the continuous curve. Note that no fitting parameter is involved. 



The scaling function ip is therefore Gaussian, 

ip(x) oc exp(-(i/+ 1 - 26)x 2 ) . (19) 

Such a solution is meaningful only if v + 1 — 26 > 0. This scaling behavior, hereafter referred to as type a (not to 
be confused with the restitution coefficient ao) is compared in Figure [7] with Monte Carlo results. The agreement 
between the analytical prediction and the numerical data involves two aspects : first, the exponent io = 1/2 allows 
us to rescale all distributions onto a single master curve. Second, this curve is exactly of the form (|19[) . where the 
prcfactor hidden in the proportionality sign cx follows from normalization. Note that the excellent agreement between 
numerical data and analytical prediction is therefore obtained without any fitting parameter. 

Type (5-scaling: On the line v + 1 — 26 = 0, the term of order q 2uJ in (|17| vanishes. If < v < 1, the terms q and 
^(i^+2)w can k e k a i ancec i w j^ n the result u = l/{v + 2), and one finds at large \x\, 

V>(z)ocex P ('--J^K+ 2 y (20) 
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FIG. 8: Scaling of type f3. Same plot as Figure [7] displaying e w /(c) vs. x = (|c| — c*)/e w at (#, f ) = (3/4,1/2) with 
= l/(f + 2) = 2/5 and c* = l/v2. The continuous curve shows the prediction of Eq. (|20p where the prefactor is determined 
from the constraint J dap = 1- 

Unlikc in a-scaling , it is not possible here to obtain ^ in close form. We will refer to (j20|) together with ui = l/(v + 2) 
as a scaling of type f3. Figure [5] shows that this behavior is in good agreement with the simulation results. Note that 
Eq. (|2T))) seems to hold also for small values of x, whereas it is a priori only valid to describe the large— a; tail of ip. 
We also note that the large— a; tail (j20f exhibits the same exponent v + 2 as in White Noise driving (see Eq. (fTTj) ). 
The important qualitative difference between the velocity distributions reported in section HT1 and here is: with WN 
driving /(c) is regular when ckq — ► 1 _ , and with Negative Friction it develops singularities. 

On the line v + 1 — 29 — 0, now with v > 1, one can only balance the terms in q and q 3uJ in (fTT)) . This leads to 
lu = 1/3, but the associated function ip diverges for large arguments, and is thus unphysical. This is an indication 
that the double peak limiting form cannot be valid on the line u + 1 — 2# = if ^ > 1. Monte Carlo results confirm 
this. They display a limiting form with three peaks in this region of parameter space (see Figure 2] and section UlIBI 
for a more thorough discussion, in particular Figure 0- 

The special point on that line with (9, v) = (1,1) represents free cooling of inelastic hard rods, which has been 
studied in [17j. There it was shown that p6[) holds, with ui — 1/3 and an asymmetric scaling function, behaving at 
large arguments as, 

ip(x) ~ Cexp [±£ 3 + o(l)] (x — > -co) 

ijj(x) ~ C'exp [-2a;( 2 / 2 ) v , +o(l)] with C = C exp [-^(y 3 )^] (as -> +oo). (21) 

This asymmetry looks quite singular since most other scaling functions identified so far are symmetric. However, 
more exceptional cases with asymmetric scaling forms tf> can be found in Figure [10] (which corresponds to the point 
(0,-1), a case of marginal stability where v = 6 — 1), as well as in the scaling shapes of the "Zoo" region of Figure 
[SJ We also emphasize that the Dirac peaks appearing at the level of description when ceo — * 1~ are not artifacts of 
discarding any spatial dependence in ([T]) , but provide the exact solution of the Boltzmann equation where due account 
is taken of the spatial degree of freedom of the particles. This has been confirmed in [j| and where the velocity 
distributions of the homogeneous Boltzmann equation has been compared with its exact counterpart, obtained by 
Molecular Dynamics simulations. 

2. Case v = (a! -scaling). 

Type a' -scaling: When v vanishes (Maxwell models), the first three terms on the rhs of (|17[) are of the same order, 
so that uj = 1/2 and ip satisfies, 

tp'(x) + 2(1 - 29)xip(x) + 2ij){x){x - = ip'(x) + 4(1 - 9)xip(x) = 0, (22) 

since (y)^ = 0. Its solution is, 

ip(x) tx exp[-2(l - ff)x 2 ]. (23) 

This scaling, coined a', a priori holds for < 8 < 1, as follows from v = and the stability requirement v+1 — > 0. 
However, when 9 = 1 (free cooling), (|23p becomes unphysical. This is a consequence of the peculiar behavior of the 
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freely cooling one-dimensional Maxwell model: the velocity distribution, which is algebraic, does not depend on ceo 
for a>o < 1 [2fj|,l2l|. In the (trivial) quasi-elastic limit, / can consequently not develop a singularity of any kind. The 
scaling ansatz p^|) has to break down for 9 = 1, and it does. For 9 < 1, the simulation data are in good agreement 
with a'-scaling predictions (see Figure [5]), where again no fitting parameter has been used. 



3. Case v < (/3-scalvng). 

Type [3-scaling: Finally we investigate the region v < 0, which is further confined by the stability requirements for 
the steady state, (9 > 0, v > 9 — 1). We then have (2 + v)oj < 2uj and the terms of order q and q^ 2+uS > UJ balance each 
other in Eq. (17]). This implies u = l/(v + 2) and 

^ = -2{\x-y\ v (x-y)U. (24) 

We recover the /3-scaling, and in particular the large— \x\ expression, 

^(x)cxexpf 2 \ X \»+A. (25) 



In the region v < 0, we have successfully tested the validity of /3-scaling against Monte Carlo simulations. In some cases 
however, we note that the best rescaling with respect to inelasticity is obtained with an exponent w opt that slightly 
differs from the predicted 1/(^ + 2). For instance, with (9, v) = (0.7,-0.2) we find w op t — 0.58 while 1/(^ + 2) ~ 0.55. 
This could indicate that the scaling limit has not yet been reached, or it reflects the fact that for negative values of 
v, it is more difficult to reach the steady state in Monte Carlo simulations. Here the collision frequency is dominated 
by encounters with \c\ — c%\ -C 1, that lead to a negligible change in the velocities of colliding partners. Consequently 
the numerical efficiency of our algorithm drops significantly. 



B. Range of validity of scaling predictions: towards a phase diagram. 

We have reported above a good agreement between DSMC calculations and the scaling predictions assuming the 
limiting double peak forms of a-,/3- and a'-scaling, when e — 1 — ao — ^ + , for several points in the {9, ^)-plane. 
We have also shown that in some (complementary) regions of this plane the scaling hypothesis (|16p does not provide 
an attracting fixed point solution of the stationary nonlinear Boltzmann equation (fT]) in the quasi-elastic limit. 
In fact, our analysis shows that physical solutions with a-scaling do not exist for v > 0, and v + 1 < 29 (e.g. 
(0,i/) = (1.1,0.3); (1.1, 1.1)), and likewise for /3-scaling with v > 1 and v + 1 = 29 (e.g. (9,v) = (1.1, 1.2)). In these 
regions we have no predictions. 

Furthermore, in the triangular region [9 > 1, v > 1, v + 1 > 29] the NF driven kinetic equation admits - at least 
on the basis of the criteria developed - QE- limiting solutions with two peaks, consistent at least to second order in s. 
The dynamics selects a different solution with three peaks (see Figure |4] with (9, v) — (1.1, 1.3)). 
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A systematic Monte Carlo investigation of various points (9, v) -shown as dots in Figure [6]- reveals that the range 
of validity of a-scaling is in fact limited to 9 < 1 |24| . In principle it would be possible to repeat the analysis of 
section Till Al with a three-peak solution, however we did not try to carry out this analysis. The reason is three-fold: 
First, it would not explain why a-scaling is not selected by the dynamics in the angular region above. Second, it 
is cumbersome. Third, there is numerical evidence that /(c) may differ from the two-peak or three-peak form (see 
the Zoo region of Figure [6]). Thus, such analysis would in any case not provide a complete picture. A numerical 
investigation appears unavoidable and was used to identify the different regions of the "phase diagram" shown in 
Figure [5] We summarize our main findings: 

• a-Scaling holds for v > and v > 29 — 1 as found analytically, with the restriction 6 < 1 that follows from 
numerical evidence. 

• /3-Scaling applies to the line < is = 29 — 1 < 1 and also to the region v < 0, where the additional restriction 
v > 9 — 1 follows from the stability requirement of the steady-state solution of Eq. fl} . Figure [TO] with 
(9, v) = (0, —1) shows a case of marginal stability where v = 6 — 1. It is observed that the scaling exponent is 
still given by w — l/(u + 2) but the form ([25]) breaks down (ip becomes asymmetric). 

• a'-Scaling is valid on the "Maxwell" line v = for 9 < 1. 

• Hard rods under free cooling, i.e. (9, v) = (1, 1), display specific scaling (see Eq. (|2ip. Tableland Ref. [lTj for 
details). 

• None of the above scalings hold for 9 > 1, where we have always observed a triple peak as in Figure [4j Figure 
1111 (with the same parameters as Figure 0]) shows that the distributions are also self-similar, with an exponent 
lu = 1/2. Although, we have no prediction for the three-peak forms, we note that the scaling function in Figure 
[IT] is compatible with Eq. (fT9|) , pertaining to a-scaling. This could be specific to the parameters chosen, since 
those value of 9 and v obey the inequalities, v > and v > 29—1, where a-scaling provides a two-peak solution 
of the Boltzmann equation. Two- and three-peak shapes therefore seem to have common features. We did not 
explore self-similarity further in the three-peak region. 

• There exists another triangular region in the (9, i/)-plane (the Zoo in Figure [6]), where /(c) does not evolve 
toward a two-peak or a three-peak form. In some instances, the Monte Carlo data are compatible with a four- 
peak limit as e — > + (see Figure [T2l where only the sector with c > has been shown). In some other cases, 
we observe precursors of what seems to be a six-peak form (see Figure [T5)) . For both Figures [T2l and [T51 it is 
difficult to decide if the limiting form for e — ► + will be a collection of Dirac distributions (with therefore a 
support of vanishing measure), or a distribution with finite support. We could however identify some points in 
the triangle where the limiting /(c) clearly is of finite support (see Figure [5]). 

To summarize, a'- and /3-scaling apply in the whole domain where they provide a solution to the Boltzmann equation, 
but a-scaling has a restricted domain of relevance compared to the region where the corresponding solution is self- 
consistent. The key features of the analytical predictions are recalled in Table UJ 



Scaling type 


rescaling exponent 


rescaling function 


a 


LU = 1/2 


ip(x) oc exp(-(y + 1 - 29)x 2 ) 


a' 


LU = 1/2 


%p[x) oc exp(-2(l - 6)x 2 ) 




w= 1/(^ + 2) 


ip(x) ~ exp(-2|xf' +2 /(^ + 2)) at large \x\ 


e = i,u = i 


LU = 1/3 


Eq. ED, ij) asymmetric 



TABLE I: Theoretical predictions for the different scaling behaviors in the two-peak region. 



IV. CONCLUSION 

We have studied the one-dimensional nonlinear Boltzmann equation in the limit of quasi-elastic collisions for a class 
of dissipative fluids where material properties are encoded in an exponent v such that the collision frequency between 
two particles with velocities c\ and ci scales like \c\ — c^Y ■ Two driving mechanisms have been considered: stochastic 
white noise (in which case the generic effects only depend on v) and deterministic negative friction (in which case, 
in addition to v, a second important exponent 9 > 0, characterizing the driving, has been considered). In both 
stochastic and deterministic cases, the quasi-elastic limit does not commute with the limit of large velocities. This is 
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C 

FIG. 10: Velocity distribution at (9, v) — (0, —1) for two inelasticities, /(c) develops a cusp at c = c* where c* is e-dependent. 
Plotting ip{x) — e/(c) vs. x = (c — c*)/e (inset), shows that here the scaling form (|16p applies with an asymmetric tp, and 
u = 1 = l/(f + 2). For such low values of the CPU time required to gather statistical knowledge is particularly large. 




-6 -4 -2 2 4 6 

x 

FIG. 11: Plots of ip(x) — e w /(c) vs. x = (c — c*)/e", with w = 1/2 for the same parameters as in Figure [4] (8, v) — (1.1, 1.3), 
i.e. for /(c) approaching a structure composed by three Dirac peaks. Here, c* ~ 0.76 denotes the position of the right peak 



seen in Figure [4] The continuous curve shows expression (|19|l . 



specific to one space dimension. There are however important differences between the two driving mechanisms. In 
the white noise case, the normalized velocity distribution /(c), suitably rescaled to have fixed variance, is regular and 
displays stretched exponential QE-tails of the form exp[— c !/+2 ]. On the other hand, f(c) with deterministic driving 
-which encompasses the much studied homogeneous cooling regime- develops singularities as e = 1 — ao — > + . 
The corresponding scaling behavior is particularly rich. We have classified the scaling forms encountered in several 
families, see Figure [5] for a global picture. Some regions of this (9, ^)-diagram are well understood, such as regions 
a, a' and /3. Some other domains resist theoretical understanding. Even if some progress might be possible in the 
three-peak region (one at c = 0, the two others at ±c*), the situation in the central Zoo region of Figure [H] seems 
more difficult to rationalize, and computationally elusive, since one needs to reach extremely small values of e to see 
the precursors of presumed singularities. 



APPENDIX A 



In this appendix, we expand the collision operator I(c\f) defined in in powers of q = (1 — oe)/2, up to second 
order. Such an expansion is required to unveil the internal structure of the singular peaks that develop as q — > + 
with driving by negative friction. Assuming that the functional form of the velocity distribution is given by (|16[) , our 
goal is to obtain here the differential equation fulfilled by the scaling function ij). Starting from 

I(c\f)= f du\c-u\ v f{u) 



(c-u) -/( C ) 



(Al) 
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FIG. 12: Two Zoo members at (6,v) = (0.75,0.3) (left) and (0,v) = (0.85,0.5) (right). Note the x- and y-scale in the plot on 
the right. For these parameters one needs to decrease e below 10 -6 to realize that a seemingly single peak is likely to split into 
two sub-peaks as e — > + . 




one obtains by extending ((9]) to 0(q 2 ) included, 



f(c) / duf(u)\c-u\ u (c-u) 



Mf)=4r-<j j ' 2 y dc 

Inserting the ansatz (fTB]) into (|A2[) . the term T^ 1 ) of order q reads 



6 \ 2 /2\ y+1 d 



4£ 
6 



b J dx 
2EY +1 d 



1 \AE J \ b J 
and expanding (1 — E(x + y)) u+1 further yields, 

v(y + l) 



/(c) / du\c-u\^f(u) 



1>(x) J dy^{y){l-E{x + y)Y +1 
ip(x) / dytp(y)\x-y\ v (x-y) 



0(q 3 ). 



(A2) 



(A3) 



7 (i) 



9 /'2V" 1 d 



-l + a;.E(i/ + l) 



AE 2 \b J dx 
The second order term 1^ in (|A2|) is 



E 2 / dyi>(y)(x + y) 2 + E» +1 / dy^(y)\x - y\ v (x - y) 



ip(x) 
(A4) 



2 V 4E J 2E \dx 



i/j(x) / dyip(y)\x ~ y\ 



u+2 



v+2 



i>(x) / dy^(y)\l - E(x + y)\ 



v+2 



(A5) 
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and we can keep only the largest term as q — > + , i.e. 



(A6) 



Collecting terms yields finally. 



q iiy- 1 d 



dx 



AE 2 \ b 



-l + xE(v + l) - 



-E 



v+l 



dyi>{y)\x - y\ v {x - y) 



(A7) 



The next step is to evaluate the right hand side of equation (|14j) . by an expansion of the moments n v and fie+i, 

l-2Ey' 



He +1 = J dc\c\ e+1 f(c) = J dy^(y) 
= b- 1 - e [l + 26{9 + l)E 2 (y 2 ) 4 ,] 



/.-,. = / / dcdu\c - u\ u+2 f(c)f(u) 

2 u+l 



b u+2 
2 u+l 

b v + 2 



dydzip{y)t/)(z)(l - E(y + z)) v+2 + E u+2 J dydzip(y)ip(z)\y - z\ u+2 
1 + {u + l){v + 2)E 2 (y 2 )^, + E»+ 2 f dydz^{y)^{z)\y - z\»+ 2 



With c close to —1/6, one can also perform the expansion 



i W" 1 /) = ~^ mil 29Ex + 29(9 - l)E 2 x% } 



(A8) 



(A9) 



(A10) 



(All) 



and the rhs of equation p^|) can thus finally be written as 



pq 



T± W^f) - ~& (ir 1 ± - ^Ex + 29(9 l)E 2 x 2 



+ ((u + l)(v + 2) - 29(9 + l))E 2 (y 2 )^ + E»+ 2 J J dydz\y - z\»+ 2 ^(y)i>(z)]} 
We obtain ifTT]). after simplification by the prefactors and integration with respect to x. 



(A12) 
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